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Steady flows that optimize heat transport are obtained for two-dimensional Rayleigh- 
Benard convection with no-slip horizontal walls for a variety of Prandtl numbers Pr and 
Rayleigh number up to Ra ~ 10®. Power law scalings of Nu ^ Ra^ are observed with 
7 « 0.31, where the Nusselt number Nu is a non-dimensional measure of the vertical heat 
transport. Any dependence of the scaling exponent on Pr is found to be extremely weak. 
On the other hand, the presence of two local maxima of Nu with different horizontal 
wavenumbers at the same Ra leads to the emergence of two different flow structures as 
candidates for optimizing the heat transport. For Pr < 7, optimal transport is achieved 
at the smaller maximal wavenumber. In these fluids, the optimal structure is a plume of 
warm rising fluid which spawns left/right horizontal arms near the top of the channel, 
leading to downdrafts adjacent to the central updraft. For Pr > 7 at high-enough Ra, 
the optimal structure is a single updraft absent significant horizontal structure, and 
characterized by the larger maximal wavenumber. 


1. Introduction 

Thermal convection is a pervasive phenomenon found in a wide variety of natural and 
engineering processes (Lappa, 2009). Examples include motion of the Earth’s molten 
core, atmospheric and oceanic dynamics and solar phenomena. A common approxima¬ 
tion in fundamental studies of convection is that all fluid parameters are constant except 
for the fluid density in the buoyancy force (Lappa, 2009; Chandrasekhar, 2013; Drazin 
and Reid, 2004) . This is the Oberbeck-Boussinesq approximation, valid in many physical 
situations, that allows systems to be described by the incompressible Navier-Stokes equa¬ 
tions subject to a buoyancy force and an advection-diffusion equation for the temperature 
held. The canonical problem is Rayleigh-Benard convection which focuses on flow fields 
driven by buoyancy and confined between two infinite horizontal plates. Rayleigh-Benard 
convection was an early study of hydrodynamic instability in which a given flow configu¬ 
ration is unstable at certain wavenumbers as a control parameter is varied. In this case, 
the control parameter is the Rayleigh number Ra characterizing the relative strength 
of buoyancy driven inertial forces to viscous forces. In atmospheric and solar environ¬ 
ments, estimates of Ra are typically exceptionally large. Caution must be exercised in 
such parameter regimes, however, because the Oberbeck-Boussinesq approximation may 
no longer be valid. 

A basic question in thermal convection is how much heat is transported for a given 
temperature difference. For the Rayleigh-Benard problem this is asking how the non- 
dimensional heat flux Nu depends on the Rayleigh number Ra. For sufficiently large Ra, 
one expects a power law scaling Nu ^ Ra'^ characteristic of a statistically self-similar 
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asymptotic regime. Transitions between power laws would indicate significant changes 
in the structure of the flow fields transporting the heat. Rigorous mathematical upper 
bounds on heat transport provide support for power law scaling Nu Ra'^ and establish 
that 7 < 1/2 for no-slip boundary conditions. Next, we briefly review the theoretical 
results emanating from upper bound theories and scaling arguments. We then briefly 
describe experimental and computational results. Comprehensive details on the history 
of Rayleigh-Benard convection can be found in various reviews (Siggia, 1994; Ahlers, 
Grossmann, and Lohse, 2009; Chilla and Schumacher, 2012). 

The original idea and the inspiration for determining upper bounds on heat transport 
is due to Malkus (1954) who sought to determine maximum heat transport subject to in¬ 
tegral constraints together with marginal stability assumptions on the mean temperature 
profile and the smallest scale of motion. Howard (1963) removed the marginal stability 
assumptions from Malkus’ theory to derive mathematically rigorous upper bounds on 
heat transport. Howard attempted to maximize Nu subject to ‘power integrals’ derived 
from the Oberbeck-Boussinesq equations and the no-slip boundary conditions. In such 
upper bound theories, the class of possible optimizing fields is expanded; that is, the solu¬ 
tions to the variational problem are not required to satisfy the full Oberbeck-Boussinesq 
equations. Howard performed his analysis with and without the continuity constraint. 
Without the constraint of continuity he found Nu^ ~ (3i?a/64) ' for large i?a, where 
Nu^ denotes the upper bound on Nu. The incorporation of the continuity equation is 
more mathematically challenging. Assuming that the optimum contained only one hori¬ 
zontal wavenumber, Howard derived Nu^ ~ . However, Busse (1969) showed 

that multi-wavenumber, divergence-free fields have higher heat transport and he derived 
1 /2 

Nug ^ (i?a/1035) ' using an approximate nested boundary layers approach. 

In the 1990s, an alternative bounding approach called the background field method 
was developed (Doering and Constantin, 1992; Doering and Constantin, 1994; Constantin 
and Doering, 1995; Doering and Constantin, 1996). The basic idea is to decompose the 
velocity field into a background field and fluctuations about the background. The back¬ 
ground fields are used as trial functions in the variational statement to maximize desired 
flow quantities (such as Nu). When applied to the Rayleigh-Benard problem with no-slip 
boundary conditions, the method yields an upper bound of Nu^^ ^ 0.167i?a^/^ — 1. Ker- 
swell (2001) demonstrated a connection between the Malkus-Howard-Busse method and 
the Doering-Constantin background field method, illustrating that the two methods ap¬ 
proach the same bound, with approach from below (above) in the Malkus-Howard Busse 
(Doering-Constantin) approach. The previous upper bounds were improved, and shown 
to be “sandwiched” between 0.031i?a^/^ (derived from the Malkus-Howard-Busse the¬ 
ory by Busse’s approximate nested boundary layers approach ) and 0.0335i?a^/^ (from 
the background field method). Although we focus on the no-slip case in this work, it 
is important to note that an upper bound for the two-dimensional free-slip case of 
7 = 5/12 < 1/2 was recently found (Whitehead and Doering, 2011). That result is 
significant because it indicates that boundary conditions and dimensionality play a role 
in limiting heat transport. Recently Wen et al. (2015) sharpened the free-slip bound and 
derived Nug — 1 < 0.106i?a®/^^. Extending the results on rigorous upper bounds to in¬ 
clude Prandtl number {Pr = i//k) effects has been challenging although some progress 
has been made in the infinite Pr limit (Whitehead and Wittenberg, 2014). 

The works discussed so far focused on developing rigorous upper bounds on heat 
transport. There are also scaling results based on classic turbulence modeling. Kraich- 
nan (1962) used eddy viscosity and mixing length ideas to develop a scaling theory of 
Nu as a function of Ra and Pr. His work led to a prediction of an ultimate regime 
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with Nu ~ and Pr-dependence as follows: Nu ~ Pr^/^Pa^/^ [In for 

Pr < 0.15; Nu ~ Pr“^/^Pa^/^ [In (Po)]~^^^ for 0.15 < Pr < 1. Later, Shraiman and 
Siggia (1990) found Nu « 0.27Pr~^^'^Ra'^^'^ assuming the existence of turbulent bound¬ 
ary layers. The more recent Grossmann-Lohse (GL) theory decomposes dissipation rates 
into bulk and boundary layer contributions (Grossmann and Lohse, 2000; Grossmann and 
Lohse, 2001; Grossmann and Lohse, 2002; Grossmann and Lohse, 2004). The Grossmann- 
Lohse framework results in two equations for Nu (Ra, Pr) and Re (Ra, Pr) with five pa¬ 
rameters that are to be fit using experimental data. The theory has also been adapted 
to the case of turbulent boundary shear layers (Grossmann and Lohse, 2011) yielding 
Nu ~ [PrRa)^^'^ x logarithmic corrections. Recent updates to the parameters based on 
experiments with aspect ratio T less than or equal to one (Stevens et ah, 2013) have 
led to excellent agreement of this theory with existing experimental data (Niemela and 
Sreenivasan, 2006; He et ah, 2012a; He et ah, 2012b). A major contribution of the GL 
theory has been the conceptual breakdown of the scaling laws into the Pr — Ra phase 
space. Our present work corresponds to regions and //„ of this phase diagram as well 
as the beginning of region IVu- In each of these regions the thermal boundary layer is 
nested within the viscous boundary layer. 

A large number of experiments have been conducted to study Rayleigh-Benard con¬ 
vection. It is not expected that the rigorous theoretical upper bounds described in the 
beginning of this section will be realized, although some experimental evidence does indi¬ 
cate that the ultimate regime may be attained (Ghavanne et ah, 2001; He et ah, 2012b). 
A main focus of experiments is on determining just how turbulent heat transport behaves 
in nature, as well as on exploring the limits of the Oberbeck-Boussinesq approximation in 
such systems (Niemela and Sreenivasan, 2006). A variety of fluids have been used as work¬ 
ing fluids in experiments (see Ghilla and Schumacher, 2012, for an overview). Rayleigh 
numbers as high as Ra « 10^^ have been reached using cryogenic helium gas (Niemela 
et ah, 2000). In most of these experiments it appears that the scaling of Nu is steeper 
than Ro?/'^ and slightly below Ra}^^ (with notable exceptions including Ghavanne et ah, 
2001, Roche et ah, 2001 and He et al. (2012b)). Indeed, Niemela et al. (2000) found 
Nu ~ We also note that the Pr dependence of these scaling laws is found to be 

very weak for Pr > 1 (Ahlers, Grossmann, and Lohse, 2009). Moreover, recent experi¬ 
mental results have been interpreted as a transition to the ultimate regime of thermal 
convection (He et ah, 2012a; He et ah, 2012b). Those results suggest that the classical 
scaling of Nu ~ Ra}^^ persists up to Ra « 10^^ after which a transition to the ultimate 
state is obtained. 

Gomputational studies of turbulent Rayleigh-Benard convection have also explored 
scaling laws for vertical heat transport (Verzicco and Camussi, 2003; Amati et al., 2005; 
Stringano and Verzicco, 2006; Verzicco and Sreenivasan, 2008). Amati et al. (2005) per¬ 
formed direct numerical simulations of Rayleigh-Benard convection within the Oberbeck- 
Boussinesq approximation up to Ra = lO^'^ at Pr = 0.7. They found a scaling of 
Nu ~ Ra}/^. Numerical studies also seek to explain large scale flow structures (Stringano 
and Verzicco, 2006) as well as discrepancies between experimental and numerical re¬ 
sults (Verzicco and Sreenivasan, 2008). The scaling laws from the computational studies 
just mentioned are derived from fully turbulent flow fields which permits contact with 
the original work by Malkus (1954) as well as with the GL theory. Several data sets from 
numerical and experimental studies can be found in Stevens et al. (2013) along with 
comparison to the GL theory. Other computational studies have explored the effect of 
Pr on heat transport. Verzicco and Gamussi (1999) and Breuer et al. (2004) performed 
simulations of Rayleigh-Benard convection in order to discern Pr effects on the flow field. 
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Verzicco and Camussi (1999) performed simulations in a cylindrical cell with diameter to 
height ratio D/H = T = 1 whereas Brener et al. (2004) performed simulations in a box 
with r = 2. Both studies found a threshold in Pr at Pr « 0.3 above which Nu displays 
very little variation with Pr. Below that threshold the flow is dominated by a large scale 
wind and Nu is highly dependent upon Pr. 

In this work, we numerically compute steady solutions to the Oberbeck-Boussinesq 
equations that optimize vertical heat transport. Thus we make contact with the upper 
bound theories mentioned above. However, our solutions exactly satisfy the Oberbeck- 
Boussinesq equations and do not include any approximations other than those introduced 
for numerical discretization. Previous work has explored such solutions at Pr = 1 (Wal- 
effe, Boonkasame, and Smith, 2015) as well as for porous medium convection (Wen, 
Corson, and Chini, 2015). Here we consider the two-dimensional problem with no-slip 
boundary conditions for Pr = 1,4,7,10,100 and Ra < 10®. We discuss the structures 
that contribute to optimal vertical heat transport. In particular, we explore the effect of 
Pr on the optimal solutions and its impact on scaling laws. A large body of early theoret¬ 
ical work considered Rayleigh-Benard convection at Pr = 7 in a domain of length 2^ ja 
where a = 1.5585 is the wavenumber to which the conduction state is unstable for no slip 
boundary conditions (Drazin and Reid, 2004). We refer to solutions with that wavenum¬ 
ber a = 1.5585 as “primary” solutions; these are the first nonlinear steady solutions that 
bifurcate from the conduction state. However the optimum heat transport solution corre¬ 
sponds to a wavenumber that generally increases with i?a, starting from the fundamental 
wavenumber a = 1.5585 at Ra = 1708. For some Prandtl numbers, we find the possible 
emergence of a 2nd local optimum transport solution and possible abrupt transitions 
in the wavenumber of the global optimum steady solution, from one local optimum to 
another. In spite of these significant jumps in optimum wavenumber for some Pr, we do 
not find any significant Pr dependence in the Nu {Ra) scaling. However, a significant Pr 
dependence is observed in determining precisely which structures lead to the optimum in 
Nu; in particular, Pr controls the length scales of these structures. Interestingly, Pr « 7 
(water) represents a demarcation between fluids that require smaller wavenumbers and 
those that require larger wavenumbers to optimize vertical heat transport. 

The remainder of the paper is as follows. In section 2 we introduce the governing 
equations, boundary conditions, and dimensionless parameters. Section 3 discusses the 
computational method used to find fixed points and optimal solutions. The main results 
are presented in section 4, followed by concluding remarks in section 5. 


2. Background 


2.1. Governing Equations 

We consider a fluid under the influence of gravity and confined between two infinite hor¬ 
izontal planes situated at y = zth so that the channel height is H = 2h. The Oberbeck- 
Boussinesq approximation is employed wherein compressibility effects (i.e. density vari¬ 
ations) are only significant in the buoyancy term. The evolution of the velocity field 
u (x, t) = {u, v) in a rectangular domain with cartesian coordinates x = {x, y) is given 

by 


^U. 

——h u • Vu -f VP = -I- a^gTy 

dt 

V • u = 0 


( 2 . 1 ) 

( 2 . 2 ) 


where P (x, t) is the kinematic pressure, v is the kinematic viscosity, is the volumetric 
thermal expansion coefficient, g is the acceleration of gravity, and y is the unit vector 
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in the vertical direction. The pressure and the horizontal velocity can be eliminated by 
taking the double curl of equation (2.1). A projection in the vertical direction then yields 
the vertical velocity equation 

The horizontal velocity is recovered from the continuity equation (2.2). The temperature, 
T (x, t), is obtained from the thermal energy equation that reduces to an advection- 
diffusion equation for the temperature, 

dT 

— + u • VT = kV^T. (2.4) 

In equation (2.4), k is the thermal diffusivity. Note that in general equations (2.3)- (2.4) 
must be supplemented by an equation for the horizontal average of the horizontal compo¬ 
nent of velocity which we refer to as the mean flow u (y, t). It is expected that the mean 
flow will decrease vertical heat transport and therefore optimal solutions for vertical heat 
transport will result when the mean flow is zero (Kerswell, 2001). We therefore directly 
set the mean flow to zero and impose mirror symmetry about the vertical plane a: = 0 
so that the solutions maintain zero-mean. Mirror symmetry is expressed mathematically 
as [m,w,T] (x,y,t) = [-u,v,T] {-x,y,t)- 

In order to fully determine the velocity and temperature helds, equations (2.3) and (2.4) 
must be augmented with boundary conditions. We impose periodic boundary conditions 
in the horizontal direction for both fields. At the top and bottom surfaces we prescribe 
uniform cold and hot temperatures, respectively, 

T {x,±h,t) = TTs, Ts>0. (2.5) 


No-slip conditions on the velocity field, m = 0, w = 0, are imposed at each surface so that 

dv 

v(x,±h,t)=0, 7^ (a:, ±h, t) = 0, (2.6) 

oy 

where the second equation follows from continuity. 

The primary diagnostic quantity is the net heat flux through the top and bottom plates 

(2.7) 

y=±h 

where the overbar (•) denotes a horizontal average. In statistically steady state, H+ = 
TL_ = TL is constant. The Nusselt number Nu is the ratio of the net heat flux TL to the 
conductive heat flux TLq = nTg/h^ thus 

_ H _ __h ^ 

KTs/h Ts dy 


( 2 . 8 ) 

=+h 


n±{t) = -f 


dT 

dy 


2.2. Nondimensionalization 

In equations (2.3)-(2.4) we nondimensionalize lengthscales by the half-height h. The tem¬ 
perature is nondimensionalized by the half-temperature difference between the bottom 
and top plates, AT/2 = Tg. The time scale is chosen as r = \fhfg' where g' = ga^Tg is 
the reduced gravity and 2r is the time that it takes a blob of fluid at temperature —Tg 
to free fall from the top wall to the bottom wall. A velocity scale is then U = h/r which 
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is the average free fall speed. The resulting nondimensional equations are 


5V 


d 


d^T 


V-u = 0 


dT 

— + u • VT = 
at 


with boundary conditions 


r(±i) = Ti 

dv 

n (a;,±l,t) = 0, — (x, ±1, t) = 0. 

ay 

This nondimensionalization leads to dimensionless diffusivities 

V , K 


= 


Vg'h^ 


and 


K* = 




(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 

(2.14) 


which are related to the Prandtl number, Pr = vj k, and the Rayleigh number, Ra = 
2g'H^/vK, through 


16Pr 


and 


K* = 


16 


RaPr 


Ra 

2.3. Heat Transport 

In our non-dimensionalization the Nusselt number, (2.8) is simply 


(2.15) 

y=-i 

Integrating (2.11) over horizontal planes, using continuity to write u • VT = V • (uT), 
yields 

assuming a statistical steady state (steady horizontal averages). This indicates that the 
net heat flux through horizontal planes is constant 



dT dT 

-K*— +vT = -K* — 
dy dy 


= K* Nu 


y=-l 


(2.17) 


since u = 0 at the bottom plate. Averaging (2.17) over the channel height, from T = 1 
at y = —1 to T = —1 at y = +1, then gives 

Nu=l + (2.18) 

Hi, 

where (•) represents a volume average. In the pure conduction state u = 0 and Nu = I. 
It is well-known, however, that the conduction state is unstable above a critical Ra of 
1708 (Jeffreys, 1928; Pellew and Southwell, 1940; Reid and Harris, 1958). Above this Ra 
the solution bifurcates from the conduction state to a steady convection state. It can 
be shown from the kinetic energy equation that (vT) > 0 whenever u 0 (again in 
statistically steady state). The onset of convection thus acts to enhance heat transport 
and Nu > 1 for Ra > 1708. In this work, we search for steady solutions that optimize 
Nu as a function of a at a given Ra. We find that Nuopt ~ RcP for sufficiently large Ra 
at fixed Pr. 














Optimal Heat Transport in Rayleigh-Benard Convection 


7 


3. Computational Method 

We seek steady, possibly unstable solutions of equations (2.9)-(2.11) that maximize 
Nu (2.18) over a. To accomplish this task, we have implemented a continuation method 
to find fixed points that uses Ra as the continuation parameter. Our implementation is 
similar to that used by Sanchez et al. (2004) and Viswanath (2007) to track periodic 
orbits in annular thermal convection and plane Couette turbulence. Net, Alonso, and 
Sanchez (2003) used a similar algorithm to find stationary states in annular Rayleigh- 
Benard convection. 

At each Ra we find steady solutions of various horizontal wavelengths, and then de¬ 
termine which wavelength results in the maximum vertical heat transport. We call this 
solution the optimal solution although it is only guaranteed to be a locally optimal steady 
solution. The richness of the Oberbeck-Boussinesq equations means that many types of 
solutions exist. For example, as will be shown in section 4.3, our optimal solutions have 
a single thermal plume emanating from the lower surface. However, it is feasible for 
other solutions to exist such as two opposing plumes with one plume emanating from 
the bottom surface and the other from the top surface. We do not expect such opposing 
plumes to result in optimal vertical heat transport but other non-steady or 3D flows 
might possibly yield higher transport. 

The length of the channel in our simulations is given by L = 2n/a where the wavenum¬ 
ber a characterizes the largest permissible wavelength in the channel. We use the flow 
map algorithm outlined in section 3.2 to find fixed points. Before discussing the fixed 
point algorithm, we briefly review the method we used for solving the fourth-order system 
(equation (2.9)). 

3.1. Solution of the Oberbeck-Boussinesq Equations 
The momentum equation for the vertical velocity field (equation (2.9)) is fourth order 
and therefore requires the four boundary conditions stated in equation (2.13). In order to 
efficiently solve these equations we follow a procedure in which an auxiliary variable, (f, 
is introduced that reduces the order of the PDF from fourth-order to second-order (Kim, 
Moin, and Moser, 1987). With (f = V^v the vertical velocity can be determined from. 


d'^T 

dx'^ 


d(j) d 
dt dx 


{u(j) — v'Tj'^u^ = -\- 


(3.1) 

(3.2) 



Although we are determining steady states, our algorithm still requires an accurate time- 
integrator (see section 3.2). In this work we use an adaptive implicit-explicit (IMEX) 
Runge-Kutta (RK) method that is third order accurate. In particular, we use a (3,4,3) 
(3 implicit stages, 4 explicit stages, third order accurate) IMEX-RK scheme (Ascher, 
Ruuth, and Wetton, 1995). Details on satisfaction of the boundary conditions and the 
time integration implementation for our particular problem are provided in appendix A. 
The spatial discretization is accomplished using a Fourier basis in the horizontal, periodic 
direction and a non-uniform finite difference scheme in the vertical direction based on 
Lagrange polynomials. The 2/3 dealiasing rule is applied in the horizontal direction to 
remove aliasing errors associated with the nonlinear terms. We checked the numerical 
resolution of our simulations in several ways. For low Ra we were able to perform a 
mesh-refinement study in which we doubled the spatial resolution. We then ensured that 
the Nu between the fine and coarse solutions agreed. In general, our simulations had at 
least 10 — 20 points in the boundary layers which satisfies the requirements suggested 
by Grdtzbach (1983). A third resolution requirement that was sometimes employed was 
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to check the horizontal energy spectra of the velocity and temperature fields to make sure 
that energy did not pile up in any wavenumbers. Finally, our results were checked against 
the Pr = 7 results of Waleffe, Boonkasame, and Smith (2015) and some of the results 
presented here for Pr = 1 and Pr = 7 were obtained with code 2 in that reference. 

3.2. Flow Map Algorithm and Fixed Points 

We now describe the flow map algorithm used to determine fixed points of this system. 
Let X = {(f), T) be a vector of solutions. We introduce a spatial discretization of the 
equations of interest (i.e. equations (3.1) and (2.11)) as described in section 3.1 from 
which we obtain a time-dependent system of ordinary differential equations, 

^=F{X). (3.3) 

One way to find steady states is to directly determine fixed points X such that 

F{X)=0. (3.4) 

Such methods represent direct-to-steady-state methods. An alternative way of deter¬ 
mining fixed points is to use the flow map. Let X be any initial state and Xt the 
corresponding solution at time T. If A is a fixed point then the difference 

Gt(A) = ^^^ (3.5) 

should be zero. Hence the goal is to find X such that 

Gt (A) = 0. (3.6) 

For a small integration time T the difference Xt — X will be small. We divide by T so 
that Gt{X) —>■ F{X) as r —>■ 0. On the other hand, it is expected that Gt (A) —>■ 0 as 
T —> oo, assuming that \Xt — X\ is bounded. Choosing T is a matter of experience for a 
particular problem, but it should be relatively small compared to the time scales of the 
modes that prevent rapid relaxation to steady state. 

We have studied the effects of the flow map approach on solving a large linear system 
of equations ^ = A A — b, where the eigenvalues of A mimic those of problems of inter¬ 
est, such as diffusion. The flow map acts as a preconditioner for the resulting Jacobian 
matrix: it can improve the condition number of the Jacobian matrix and redistribute the 
eigenvalues such that they cluster around the origin. These observations are consistent 
with previous studies (Sanchez et ah, 2004) and indicate improvements for solving some 
linear systems of equations. In particular, the generalized minimum residual (GMRES) 
method is known to converge well when the eigenvalues of the matrix are clustered around 
the origin (Trefethen and Bau III, 1997). 

3.3. Details of the Flow Map Algorithm 

Equation (3.6) represents a large system of equations with approximately N unknowns 
where N increases with Ra. In this study N « 10® for the largest values of Ra used. 
Equation (3.6) is typically solved using some variant of Newton’s method. The full ver¬ 
sion of Newton’s method can be difficult to implement and expensive to use. Difficulties 
include forming the full Jacobian matrix and efficiently solving the resulting matrix sys¬ 
tem. Sophisticated approaches exist to help make Newton’s method more cost-effective. 
These include preconditioning the Jacobian matrix, inexact Newton methods based on 
Krylov subspaces, and Jacobian-free methods (Knoll and Keyes, 2004). 

In the present work, we solve (3.6) by means of a Jacobian-free Newton-Krylov method. 
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Newton’s method results in the linear system, 

AX = - (3.7) 

X=Xk 

where Xk represents the solutions at nonlinear iteration k and AX = Xi~+i — X^. Rather 
than forming the full Jacobian, a forward finite difference scheme is used to approximate 
the action of the Jacobian on AX (Knoll and Keyes, 2004). The GMRES method is 
then used to solve the resulting linear system (Saad and Schultz, 1986; Trefethen and 
Bau III, 1997). We note that it is not necessary to form the full Jacobian (or even its 
approximation) because the GMRES algorithm only requires the action of a matrix on 
the solution. Each iteration of the GMRES method solves equation (3.7) with 


dGr 

"d^ 


dGr Gr (Xfc + eAX) — Gt (Xfc) 

^ X.X. “ ^ 


(3.8) 


with the choice 


e 


■\/^niachine 

IIAXIP 


max(Xfc • AX, ||AX||) 


(3.9) 


where Cmachine is the machine tolerance and || • |j is the L 2 norm. Several possibilities exist 
for the choice of e (see Knoll and Keyes, 2004, pg. 363) but the main point is that it 
must be small enough to get an accurate derivative in (3.8) but not so small that floating 
point errors pollute the solution (Dennis and Schnabel, 1996, section 5.4). 

We conclude this section with a few comments and observations based on our expe¬ 
rience with the flow map algorithm. The fact that the flow map algorithm redistributes 
the eigenvalues of the Jacobian matrix also has benefits when using the GMRES method 
to solve the linear system (3.7). To understand the latter, recall that the basis vectors 
for the Krylov subspace are formed from a power method, and that the power method 
focuses on the largest eigenvalue of the system. Since the flow map algorithm results in 
a larger separation between the largest eigenvalue and the other eigenvalues, the power 
method converges faster. Hence, a smaller Krylov subspace is required for the linear 
solves. 

In terms of computation time, there appears to be an optimal integration time T that 
minimizes the wall-clock time of the simulation. This optimal time is dependent upon the 
problem at hand. For example, increasing T causes all time-integrations to take longer. 
On the other hand, increasing T allows for better preconditioning and results in a much 
reduced time for solving the linear system. For the current problem, we have observed that 
T « 0.6 results in the quickest computation time for low Ra around 10®. We have used 
T = 2.0 throughout for both high and low Ra. Some of the results presented herein for 
Pr = 1 and Pr — 7 for Ra > 10® were computed with code 2 in Waleffe, Boonkasame, 
and Smith (2015) which is based on spectral Chebyshev integration and does not use 
GMRES. 


4. Results 

In this section, we present steady solutions to the Oberbeck-Boussinesq equations that 
maximize heat transport. We consider Prandtl numbers Pr ranging from 1 to 100 and 
note that these values roughly correspond to real fluids. For example, Pr = 7 is often 
used for water and engine oil has Pr of roughly 100. In subsection 4.1, we explore the 
effects of both Pr and the horizontal wavenumber a on optimal heat transport. Sub¬ 
section 4.2 investigates the impact of Pr on scaling laws, while subsection 4.3 discusses 
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a 

Figure 1: Variation of Nu with horizontal wavenumber a at Ra = 3x 10® for Pr = 1,4,10. 


the coherent structures that give rise to the optimal solutions. The mean temperature 
profiles associated with the optimal solutions are considered in subsection 4.4. 


4.1. Multiple Local Maxima and Pr Effects 

The variation of heat transport Nu with wavenumber a is shown in figure 1 for Pr = 1, 
Pr = 4 and Pr = 10 at fixed Ra = 3 x 10®. We note that two local maxima emerge 
for these three fluids. Thus, there are two scales that locally maximize heat transport. 
The solution corresponding to the global maximum, hereafter called the optimal solu¬ 
tion, occurs at different length scales depending upon the fluid under consideration. For 
example, the global maximum occurs at a small wavenumber a for Pr = 1, whereas it 
occurs at a larger ol for Pr = 10. We have indicated in the figure the ‘first’ and ‘second’ 
local maxima for Pr = 10. When we refer to the first local maximum, we have in mind 
the maximum that occurs at smaller wavenumber a, whereas when we refer to the second 
local maximum, we mean the maximum that occurs at larger wavenumber a. 

An intriguing behavior is illuminated upon inspection of other Ra and Pr. As indicated 
in figure 2, it is not always the case that two local maxima are present. Four different 
scenarios are described for Pr = 1 (2a), Pr = 7 (2b), Pr = 10 (2c) and Pr = 100 (2d). 

Figure 2a; Pr = 1; 2 x 10^ < Ra < 3 x 10®: For Ra < 10® only one maximum 
exists. However, for Ra > 2 x 10® we observe two local maxima. The 
first maximum is the optimal solution and there is no indication (at 
least for Ra < 3 x 10®) that the second local maximum will become 
the global maximum. 
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Figure 2b; Pr = 7; 5 x 10^ < Ra < 4 x 10®: For Ra < 10® there is only 
one local maximum. However, for Ra > 10® two local maxima are 
present. Moreover, the second local maximum is growing, and by 
Ra = 3 X 10®, the two local maxima are nearly identical, although 
they occur at different scales. In this i?a-range, it is not clear whether 
the first local maximum or the second local maximum will ultimately 
be the optimal solution, and the issue must be settled by computing 
maxima at higher Ra. We have performed these computations, and 
the inset presents the Nu — a behavior at Ra = 1 x 10®. We observe 
that at this Ra the second local maximum has receded compared to 
the first local maximum. 

Figure 2c; Pr = 10;5xl0‘^< Ra < 2 x 10®: Yet another scenario is observed, 
wherein the new feature is that the second local maximum actually 
becomes the global maximum above a certain Ra. The second local 
maximum appears at Ra « 9 x 10"^ and grows relative to the first 
local maximum. At Ra « 2 x 10®, the second local maximum is 
equal to the first local maximum. Beyond this i?a, the second local 
maximum becomes the global maximum. 

Figure 2d; Pr = 100; 2 x 10'^ < Ra < 5 x 10®: The large Pr = 100 case introduces 
a fourth type of behavior with Ra, where the first local maximum 
vanishes for high enough Ra. Once again, a single local maximum 
is present for Ra below a certain value, this time Ra = 8.5 x 10"^. 

A second local maximum appears at Ra = 8.5 x 10^ and begins to 
grow relative to the first local maximum. By Ra = 10®, the first and 
second local maxima are equal. Interestingly, by Ra = 2 x 10® the 
first local maximum has vanished and the second local maximum is 
the global maximum. 

For ease of discussion below, we use the term ‘coexistence interval’ to refer to the Ra- 
interval in which two local maxima are present, recognizing that the interval changes 
with Pr. 

Figure 3 summarizes the different scenarios we observed in the Prandtl number range 
1 < Pr < 100, for Rayleigh numbers Ra < 10^. The figure shows the ratio of the second 
local maximum to the first local maximum Nu^/Nu.^^ as a function of Ra. For Pr = 1, 
Pr = 4 and Pr = 7 the second local maximum does not overtake the first maximum, 
whereas for Pr = 10 and Pr = 100 the second maximum eventually becomes the global 
maximum. For Pr = 7, the two maxima are almost equal with the second local maximum 
becoming 99.6% of the first local maximum at Ra « 2 x 10® before slowing dropping off. 
We also note that the coexistence region of the two maxima depends heavily on Pr. This 
is most clearly observed for the Pr = 100 case, in which the two maxima coexist only 
for a very small interval of Ra. For Pr = 4, Pr = 7 and Pr = 10 we did not observe 
termination of the coexistence region. However, for Pr = 1 and Pr = 100 the coexistence 
region was found to be finite, with the second (first) maximum disappearing for Pr = 1 
{Pr = 100). Figure 3 implies the existence of two different parametric regions for optimal 
heat transport roughly partitioned by Pr = 7. For Pr > 7, the second local maximum 
overtakes the first local maximum to become the optimal solution and appears to remain 
the optimal solution. Since the second local maximum occurs at larger a, our calculations 
suggest that optimal solutions occur at relatively smaller scales for fluids with Pr > 7 
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Figure 2: Variation of Nu with horizontal wavenumber a for various Ra at (a) Pr = 1, 
(b) Pr = 7, (c) Pr = 10 and (d) Pr = 100. The first and last Ra contour values are 
marked in each plot. Note that the color scheme between different plots is different as is 
the Nu scale, however the a scale is identical. 


at high-enough Ra. Conversely, for Pr ^ 7, the second local maximum never overtakes 
the first local maximum, and therefore the first local maximum remains the optimal 
solution. Thus, for fluids with Pr < 7, optimal solutions occur at the larger of the two 
length scales. The Pr dependence is significant in determining the flow structure that 
leads to optimal solutions. 

4.2. Scaling Laws and Impact of Pr 

Figures 4 and 5 present log-log plots of optimal Nu and optimal a with Ra. In figure 4 we 
observe the dependence of the optimal Nu on Ra in the form Nu oc Ra^ for Ra > 10®. 
The scaling law of Nu = 0.115 RaP'^^ was determined for Pr = 7 in Waleffe, Boonkasame, 
and Smith (2015), and is shown by the dashed line in figure 4 for comparison. This scaling 
was obtained from data in the range 10^ < Ra < 10®. We observe that there is very little 
variation from this scaling with Pr although the prefactor for Pr = 1 is higher. 

Figure 5 shows the dependence of the two maximal wavenumbers ^ on the Rayleigh 
number Ra. We begin our discussion by focusing on the open markers in figure 5 (first 
maximum). For Pr = 100, the first local maximum ai vanishes after Ra = 1.3 x 10® 
and hence a clear scaling of ai with Ra cannot be determined. However, for the other 
Pr considered in this work, a scaling regime appears. The scaling of ai determined 
in Waleffe, Boonkasame, and Smith (2015) for Pr = 7 is shown in the dashed line: 














Nu2/Nu\ 
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Figure 3: Ratio of local maxima as a function of Ra for various Pr. 



Figure 4: Optimal Nu vs. Ra; global optimal data is shown and may correspond to 
distinct wavenumber. 
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Figure 5: Maximal a vs. i?a; open markers indicate the first maximum and filled markers 
indicate the second maximum; the upper dashed-dot line represents the upper bound 
from linear stability theory. 


ai = Once again, the scaling shows almost no change with Pr. For Pr = 1 

the power law of the scaling appears consistent with the other Pr although it is likely 
that the prefactor is larger. We note that for Ra < 10®, the Ra dependence of ai exhibits 
some wobbles, conjectured to be linked to the appearance/disappearance of coiling arms 
emanating from the central temperature plume of the optimal steady solutions (Waleffe, 
Boonkasame, and Smith, 2015). 

Next we turn attention to the closed markers in figure 5 corresponding to the second 
maximum with wavenumber 02 . A scaling regime emerges at Ra « 10®, and the scaling 
regime for these solutions is more complete because all fluids considered in this work 
exhibited a second local maximum in vertical heat transport that persisted for several 
decades of Ra. For Pr = 7 we find 02 ~ 0.257i?a®-^®® with the current range of Ra 
sampled by our data. It can be shown from linear stability theory {Nu = 1) that the 
upper wavenumber limit for marginal stability scales as ~ 0.5Ra^^^ for large Ra. 
This scaling is marked in the dash-dot line in figure 5. Nu drops to 1 as a approaches 
from below and there is no convective steady state for a > a„ (see also figure 2). We 
expect that data at higher Ra would reduce the best-fit 02 to be at or below 02 ~ RaS'^^. 

Figure 6 presents Nu compensated by 0.115 which is the best least square fit 

scaling of Nu at Pr = 7 observed in figure 4. In particular, we consider the Nu vs. 
Ra scaling at various Pr, and for each local maximum of Nu (first and second local 
maxima). We observe that the Nu{Ra) scaling laws that result from tracking the first or 
second local maxima are not very different. The most obvious discrepancy in Nu between 
the first and second local maxima occurs for Pr = I. However, the compensated plots 
are both flat over their respective Ra ranges. This indicates that the scaling exponent 
7 = 0.31 is observed for both types of solutions at Pr = 1. The fact that the scaling 
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Figure 6: Nu compensated by 0.115 for various Pr. Open markers: solutions from 

the first maximum. Filled markers: solutions from the second maximum. Optimal solu¬ 
tions occur when one curve is above the other curve. 


of the optimal solution at Pr = 1 is greater than unity can be attributed to a different 
prefactor in the scaling law (greater than 0.115). Similarly, the solution at the second 
local maximum has a prefactor less than 0.115. At Pr = 7 and Pr = 10 there is almost 
no discernible difference in the scaling laws between to the two types of solutions. When 
Pr = 100 the optimal solution also closely follows the scaling shown in figure 4. 

4.3. Coherent Structures and Optimal Heat Transport 

Sections 4.1 and 4.2 showed that the A^rt-maximizing wavenumbers are larger than the 
wavenumber a = 1.5585 corresponding to the onset of convection (the primary wavenum¬ 
ber). Thus the maximizing solutions have smaller horizontal scales than their non-optimal 
counterparts in the primary box. Here we compare the details of the optimal and non- 
optimal coherent structures. 

In figure 7, contours of temperature in the channel are plotted for Pr = 100 at Ra = 
1X 10^. The top plot is the non-optimal temperature in the primary box with wavenumber 
a = 1.5585. The main features are the thick horizontal temperature arms emanating from 
the main plume that roll up with increasing Ra. By contrast, the plumes in the bottom 
two plots maximize vertical heat transport and exhibit much less, if any, horizontal 
structure. In fact, the thermal plume corresponding to the optimal solution (bottom 
right) does not have any spiral arms or lobes; the heat transport is almost exclusively 
straight up and down in the vertical direction. Similar results are shown in figure 8 for 


































16 


D.Sondak, L.M. Smith and F.Waleffe 



X 




Figure 7: Temperature contours at Ra = 1 x 10® and Pr = 100. The top figure presents 
the primary solution (a = 1.5585, Nu = 4.429), the bottom left figure represents contours 
at the first local maximum (a = 3.0, Nu = 4.819), and the bottom right figure represents 
contours at the second local maximum {a = 4.3, Nu = 4.821). Note that the second local 
maximum is the optimal solution at the Ra considered. 


Pr = 10 at Ra = 2 x 10®. Again the optimal temperature structure is a single plume 
without significant horizontal structure (right; second local maximum), while the non- 
optimal plume (left; first local maximum) has coiling arms emanating from the central 
plume. Figures 7 and 8 show representative optimal solutions in fluids with Pr > 7 
at large-enough Ra, for which the optimal solution occurs at small scales (large a; the 
second maximum). 

However, as shown in figure 9, the temperature structure of optimal solutions for 
Pr < 7 seems to contradict the simple optimal structure seen for Pr > 7. The figure 9 
contours are for Pr = 1 at Ra = 5 x 10®, with optimal solution on the left (a = 3.7; 
first maximum), and non-optimal solution on the right {a = 8.1; second maximum). As 
expected, horizontal structure is prohibited by the large-a maximum, but in this case 
there is relatively weak vertical heat transport associated with the central plume. The 
smaller wavenumber (larger horizontal scale) allows for downdrafts to the left and right 
of center, which in turn allows for more vigorous updraft of heat in the center, and thus 
an overall larger vertical heat transport. Contours of the optimal vertical velocity (not 
shown here) confirm that the vertical velocity is largest in the middle of the channel. 
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Figure 8: Temperature contours at Pr = 10 and Ra = 2 x 10® corresponding to optimal 
Nu occur at a = 10.45 and Nu = 11.391 (right plot). Temperature contours for the first 
local maximum of Nu occur at a = 3.35 and Nu = 11.009 (left plot). 


Thus, depending on the Pr and Ra numbers, there is a delicate balance determining 
whether or not the optimal solution will have horizontal structure in the form of coiling 
arms emanating from the central plume. For Pr > 7 at large-enough Ra, the optimal 
solution is the ‘second maximum’ corresponding to small horizontal scale that does not 
admit significant horizontal structure. In this case, top to bottom heat transport through 
the central plume is the most efficient conduit for overall vertical heat transport. For 
^ optimal solution is the ‘first maximum’ corresponding to a larger horizontal 

scale, wherein the largest overall vertical heat transport is achieved through a vigorous 
central updraft accompanied by adjacent downdrafts. 

We close this subsection with a visualization of the optimizing solution at Ra = 10® 
for Pr = 1 and Pr — 10. (figures 10 and 11). The two leftmost plots in figure 10 present 
optimal temperature contours and corresponding streamlines at Pr = 1. The rightmost 
plots show optimal temperature contours and corresponding streamlines at Pr = 10. 
Figure 11 shows the top (top plot) and bottom (bottom plot) portions of the optimal 
temperature field with streamlines superposed at Pr = 10. A main feature to note is 
that the streamlines deeply penetrate the boundary layer. We remark that Hassanzadeh, 
Chini, and Doering (2014) found optimal solutions for a related problem in which the 
class of solntions was expanded to include all divergence free velocity fields with fixed 
enstrophy and free-slip boundary conditions. The aspect ratio in that work was found to 
scale as Ra^^^ as in the present work. One notable difference is that we do not observe 
the formation of a recirculation region near the top of the channel. 
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Figure 9: Temperature contours at Pr = 1 and Ra = 3 x 10® corresponding to optimal 
Nu occur at a = 3.7 and Nu = 7.058 (left plot). Temperature contours for the second 
local maximum of Nu (right plot) occur at a = 8.1 and Nu = 5.473. 


4.4. Mean Temperature Profiles and Nu {Pr) 

We now consider mean temperature profiles of the optimal solutions introduced in sec¬ 
tion 4.3. Figures 12 and 13 present mean temperature profiles at Pr = 10. First, in 
figure 12 we compare mean temperature profiles at Ra = 2 x 10® of the primary solution, 
the optimal solution, and the solution that gives the first maximum. The primary solution 
is nearly marginally stable, possessing a nearly isothermal region though the core of the 
channel. However, a very slight stable stratification of the primary solution can be ob¬ 
served in the inset of figure 12 just above y = 0.5. The solution corresponding to the first 
maximum is marginally stable in the center of the channel but exhibits stably stratified 
transition regions between the center and the unstably stratified boundary layers. The 
optimal solution, on the other hand, is unstably stratified in the center of the channel, 
but maintains the stably stratified transition regions and the unstable boundary layers. 
We note that non-monotonic temperature profiles have been considered in the infinite 
Pr case (Doering, Otto, and Reznikoff, 2006). Interestingly, each solution presented in 
figure 12 is close to the marginally stable configuration constructed by Malkus (1954) 
and resulting in the prediction of Nu ^ Raf^^. 

As shown in figure 4 the scaling of the optimal solution Nu ~ is less than 

Nu ^ Ra^/^ which might be attributed to small departures from marginal stability. 
Figure 13 shows optimal mean temperature profiles at Pr = 10 for 2 x 10® < Ra < 1 x 10®. 
Each of these profiles exhibits unstable stratification in the bulk and in the boundary 
layers, but still possesses the stably stratified transition regions. Perhaps surprisingly, 
the optimal solutions become more unstably stratified in the bulk as Ra increases, and 
thus departure from marginal stability seems to grow. The latter is related to the single 
updraft nature of the optimal solutions, which have increasing horizontal wavenumber 
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Figure 10: Ra = 10®. Two leftmost plots are temperature contours and streamlines for 
the optimal solution at Pr = 1 with a = 9.5 and Nu = 41.49. Two rightmost plots are 
temperature contours and streamlines for the optimal solution at Pr = 10 with a = 27.86 
and Nu = 35.63. As usual y G [—1,1] and x G [—Tr/a, 7r/a]. 


and diminishing horizontal structure; there are no significant horizontal arms or coils to 
bring colder fluid beneath warmer fluid. 

A different trend is observed for Pr < 7, as illustrated in figure 14, showing horizontally 
averaged optimal temperature profiles at Pr = 1 for 10® < Ra < 10®. At first glance, 
these profiles all appear to be marginally stable. The inset in figure 14 shows that this 
observation is not quite the entire story. Indeed, near the center of the channel the profiles 
are nearly marginally stable, with a tendency toward stable stratification, especially at 
lower Ra. Stably stratified transition regions between the center of the channel and the 
boundary layers are once again present. The mean stable stratification is linked to the 
two downdrafts on either side of the central updraft, which bring colder fluid underneath 
warmer fluid in a stable configuration. 

Figures 12, 13 and 14 suggest that for Pr <7 optimal solutions have a stably stratified 
bulk region whereas the optimal solutions for Pr > 7 have an unstably stratified bulk 
region. Figure 15 provides more evidence of this emerging picture, presenting optimal 
solutions at Ra = 5 x 10® for Pr = 1,4, 7,10,100. Indeed, the optimal mean temperature 
profiles for Pr = 1,4,7 are stably stratified near the middle of the channel, while those 
for Pr — 10 and Pr — 100 are unstably stratified. The stably stratified transition regions 
between the middle of the channel and the boundary layers are present at all Pr. We 
summarize the three regions observed in the optimal temperature profiles as follows: 

(а) The temperature boundary conditions produce unstable boundary layers. 

(б) The central region of the channel may be stably stratified (Pr ^ 7) or unstably 
stratified (Pr >7). 

(c) The transition regions between the channel center and the boundary layers were 
always observed to be stably stratified in the optimal transport solutions. 
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Figure 11: Streamlines of the optimal solution at Pr = 10 and Ra = 1 x 10® near the 
top (top plot) and bottom (bottom plot) superposed on contours of temperature where 
X G [—7r/a,7r/a] and a = 27.86. 


Finally, we compare the two dimensional optimal solutions found in the present work to 
the Grossmann-Lohse (GL) theory of turbulent convection (Grossmann and Lohse, 2000) 
and to fully three-dimensional numerical simulations (Verzicco and Gamussi, 1999). In 
figure 16 we show the variation of Nu with Pr for 0.5 < Pr < 100 at Ra = 5 x 10^. The 
work presented here corresponds to regions and //„ as well as the beginning of region 
IVu in the Pr — Ra phase diagram of the GL theory (Grossmann and Lohse, 2001). In 
each of these regions the thermal boundary layer is nested within the viscous boundary 
layer. We remark that the 2D optimal solutions are in agreement with both the GL theory 
as well as the 3D turbulent solutions. This intriguing observation begs the question: do 
the 2D steady solutions computed in the present work play a role in the dynamics of 
3D turbulent convection? For future consideration of the possible connection, one must 
consider aspect-ratio effects as well as the difference in dimensionality (2D vs. 3D). The 
theoretical Rayleigh-Benard problem has an essentially infinite aspect ratio, whereas the 
aspect-ratio dependent parameters in the GL theory have been fit to experiments with 
aspect ratio F < 1. The turbulent simulations of Verzicco and Gamussi (1999) were 
performed in a cylindrical cell with F = 1. There is evidence, however, that aspect ratio 
dependence may be weak (Ahlers, Grossmann, and Lohse, 2009). On the other hand, the 
sensitivity of Nu to the large scale circulation arising in finite aspect ratio configurations 
remains an open issue (Ahlers, Grossmann, and Lohse, 2009). 
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Figure 12: Horizontally-averaged temperature profiles for Pr = 10 at Ra = 2 x 10®. 



Figure 13: Horizontally-averaged optimal temperature profiles at Pr = 10 for 2 x 10® < 
Ra<lx 10®. 
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Figure 14: Horizontally-averaged optimal temperature profiles at Pr = 1 for 10® < Ra < 
10 ®. 
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Figure 15: Horizontally-averaged optimal temperature profiles at Ra = 5x 10® for Pr = 1, 
Pr = 4, Pr = 7, Pr = 10 and Pr = 100. 
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Figure 16: Comparison of 2D optimal solutions to the Grossmann-Lohse theory and to 3D 
turbulent solutions (Verzicco and Camussi, 1999) for 0.25 < Pr < 100 at Ra = 5 x 10®. 

5. Conclusions 

Optimal transport solutions for 2D Rayleigh-Benard convection have been computed 
for Ra up to 10® for Pr > 1. This work is an extension to that performed by Waleffe, 
Boonkasame, and Smith (2015) at Pr = 7. Most of the solutions in the present work have 
been found using a flow map algorithm in conjunction with a Jacobian-free Newton- 
Krylov solver to find steady solutions to the Oberbeck-Boussinesq equations. Some of 
the high Ra, Pr = 1 and Pr = 7 solutions were computed with code 2 in Waleffe, 
Boonkasame, and Smith (2015). An optimal solution at a given Ra and Pr represents 
the steady solution that maximizes heat transport among all horizontal wavenumbers. 
One goal is to make contact with rigorous upper bound theory (Howard, 1963; Busse, 
1969; Doering and Constantin, 1996) by computing optimal solutions constrained exactly 
by the Oberbeck-Boussinesq equations. Another intention is to compare scaling laws 
derived from experiments (Niemela et ah, 2000; He et ah, 2012a) and computational 
studies (Amati et ah, 2005) of turbulent Rayleigh-Benard convection with our optimal 
solutions in the spirit of Malkus’s original idea that turbulence maximizes heat transport 
subject to some marginal stability constraints (Malkus, 1954). A significant portion of 
this work is aimed at exploring the effect of Pr on scaling laws and optimal solutions. 

We have considered Pr = 1,4,7,10,100 and our results indicate that there is little 
to no variation of scaling of iVu with Pr (for the parameter ranges considered). This 
is consistent with studies in turbulent Rayleigh-Benard convection which report that 
Nu (Pr) is relatively constant for Pr > 1 (see Ahlers, Grossmann, and Lohse, 2009, 
figure 5). The optimal solutions exhibit a scaling of the form Nu ~ which is 

in agreement with that found by Waleffe, Boonkasame, and Smith (2015). Moreover, 
horizontally averaged profiles of the temperature field are close to a marginally stable 
profile, roughly consistent with the prediction by Malkus (1954) leading to Nu ~ Ra}^^. 
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However, the lack of Pr dependence for scaling relations is not the entire story. Indeed, 
Pr is shown to have a significant impact on the scale of the optimal solution. When 
considering the variation of Nu with a, we observe that, above a certain i?a, two local 
maxima exist; one at a smaller scale and the other at a larger scale. There are two regions 
partitioned by Pr « 7: for Pr > 7 optimal transport solutions occur at smaller scales 
and for Pr < 7 optimal transport solutions occur at larger scales. 

One might expect optimal heat transport to be realized by direct ‘wall-to-wall’ trans¬ 
port of hot and cold fluid in narrow plumes with little horizontal structure and no re¬ 
turning downdrafts of hot fluid and updrafts of cold fluid. Inspection of the optimizing 
solutions reveals that this structure is associated with the ‘second maximum’ at small 
scales, and optimizes heat transport for Pr > 7 at large-enough Ra. For Pr < 7, the 
optimal heat transport occurs at the ‘first maximum’, and is associated with a larger 
horizontal scale. In this case, the vigorous central updraft is accompanied by adjacent 
left/right downdrafts, leading to a slightly stable interior mean temperature profile. Evi¬ 
dently, this is another way to minimize the average thermal boundary layer to maximize 
transport. Alternatively, with free-slip boundary conditions and an expanded class of ve¬ 
locity fields, optimal transport solutions were found to exhibit recirculation near the top 
of the channel (Hassanzadeh, Chini, and Doering, 2014). Thus there is a delicate balance 
when optimizing over all possible velocity fields: the velocity must be strong enough to 
efficiently transport heat vertically, but not so strong that downdrafts of heat overcom¬ 
pensate for the vertical transport. It is unclear why Pr « 7 should be the demarcation 
between smaller-scale optima consisting of a central updraft only, and larger-scale optima 
admitting some horizontal structure (downdrafts adjacent to the central updraft). 

An obvious extension of the present work is to consider a wider range of Pr and Ra 
numbers and in particular Pr < 1. Stricter spatial and temporal resolutions will be neces¬ 
sary, requiring algorithmic modifications including possibly code parallelization. Another 
direction of our future work is understanding the connections between the optimal steady 
solutions and turbulent flows. It has been conjectured that such unstable solutions ac¬ 
tually control the dynamics of turbulent flows (Waleffe, 2009) and the close agreement 
between the Nu{Ra) data for optimum 2D solutions and 3D turbulent data is most 
intriguing. Is it possible to detect the signature of the optimal solutions within the tur¬ 
bulence data? One step in this direction involves performing a fully resolved simulation 
at high Ra and collecting local snapshots of the fields at times and locations where the 
Nusselt number achieves its maximal values in order to ‘educe’ a characteristic optimum 
transport structure (Stretch, 1990). These snapshots can be compared to the coherent 
structures educed through other flow decompositions (e.g., Smith, Moehlis, and Holmes, 
2005; Schmid, 2010). Thus one may begin to explore statistics such as how often the 
optimal structures are sampled by the turbulent flow. More sophisticated analyses of the 
state space and the role of the stable and unstable manifolds of the optimum solutions 
should be easier to perform than in 3D shear flows (Gibson, Halcrow, and Cvitanovic, 
2008). We hope that such investigations will be the subject of forthcoming publications. 
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Appendix A. Boundary Conditions and Time Integration 

A.l. Time Integration 

We now apply the (3,4,3) IMEX-RK method developed by Ascher, Ruuth, and Wet- 
ton (1995) to the Rayleigh-Benard problem governed by equations (3.1) and (2.11). We 
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rewrite the equations concisely as 


d ,, ^ 2 ^ d^T 

- + (j) + 2 

; ax ax^ 

(Al) 

at 

(A 2) 

= {(p, u, v) = up — vW'^u 

(A3) 

= Afr {T, u,v) = u ■ VT. 

(A 4) 


where 


We will treat the diffusion terms implicitly and all other terms explicitly. We denote by 
(/)" the field at the current time step. To evolve (p to step n + 1 we have, 




where 


t=i 




At stage f 0 is given by 

(1 - a"Ati/*V2) 0* = (^" + At 
The implicit operator is given by, 


i=i 


11 = I- a"Ati/*VA 


(A 5) 

(A 6) 
(A 7) 
(A 8) 

(A 9) 

(A 10) 


The right hand side of equation (A 9) represents in equation (A 17a) in appendix A.2. 
The temperature at step n + 1 is given by, 


cs 

T”+i = T” + At ^ b> {k{, + 


i=i 


and 


4 - 


Ktr = i ^ 1 


ifi = -A7p-i. 

The temperature at stage i is solved from the implicit equation 


(1 - a“AtK*v2) = T^ + At 


E 

i=i 








(All) 


(A 12) 
(A 13) 
(A 14) 


. (A 15) 


The values of a*l, and F can be found in Ascher, Ruuth, and Wetton (1995). 
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A.2. Boundary Conditions 


In order to satisfy the boundary conditions at stage i we let 

v’' = Vp-\- Civl + C2V2 


(A 16) 


and solve three separate problems for u®, vl, and V 2 at each stage of the time integration 
where is the particular solution and Vi and V 2 are two homogeneous solutions. The 
three problems are: 



(Al7o) 


v\; = 4, ^;;(±i) = o. 
^;i(±l) = 0. 
V^vl2=4>2, ^^2(±i) = o. 


(A17&) 


(A 17c) 


is the implicit time integration operator at stage i and contains contributions from 
the explicit portion of the time integration at stage i (see appendix A.l). The constants 
Cl and C 2 are determined so that equation (2.13) for the boundary conditions is satisfied. 
At each wall we have 


This results in a 2 x 2 matrix system that can be solved for ci and C 2 , 


dvl dvn 



r 9vi j 



Cl 



dv] , , dvn , , 




dvl, 



C2 




In the non-adaptive version of our code, v{, ^1 and (p^ are computed as part of a 

preprocessing step and used in the determination of the constants. When using adaptive 
time-step sizes v\, v\, p\ and p\ computed just prior to the next time-integration. 








